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Abstract 

We investigate the site percolation transition in two strongly correlated systems in three dimen- 
sions: the massless harmonic crystal and the voter model. In the first case we start with a Gibbs 
measure for the potential, U = ^ Z]<a;j/> — </)(y))^, x,y £ 1? , J > and G M, a scalar 

height variable, and define occupation variables ph{x) = 1, (0) for (j){x) > h{< h). The probability 
p of a site being occupied, is then a function of h. In the voter model we consider the stationary 
measure, in which each site is either occupied or empty, with probability p. In both cases the trun- 
cated pair correlation of the occupation variables, G{x — y), decays asymptotically like \x — y\~^. 
Using some novel Monte Carlo simulation methods and finite size scaling we find accurate values 
of pc as well as the critical exponents for these systems. The latter are different from that of 
independent percolation in d = 3, as expected from the work of Weinrib and Halperin [WH] for 
the percolation transition of systems with G{r) ~ r^" [A. Weinrib and B. Halperin, Phys. Rev. B 
27, 413 (1983)]. In particular the correlation length exponent v is very close to the predicted value 
of 2 supporting the conjecture by WH that u = ^ is exact. 

PACS numbers: 05.70.Jk,05.50.+q,64.60.Ak,64.60.Fr,87.53.Wz,89.75.Da 



*Electronic address: marinov@physics.rutgers.edu 
'''Electronic address: lebowitz@math.rutgers.edu 



1 



I. INTRODUCTION 



A translation invariant ergodic system of point particles on a lattice, say Z*^, in which 
each site is occupied with probability p, < p < 1, is said to percolate when it contains an 
infinite cluster of occupied sites, connected by nearest neighbor bonds. This event satisfies 
the zero/one law, i.e. the probability that the system percolates is either zero or one [1, 2]. 
For the case in which the sites are independent the transition from the non-percolating state 
for p < Pc and the percolating one for p > pc is one of the simplest examples of critical 
phenomena. The probability that a given site, say the origin, is connected to infinity, i.e. is 
part of the infinite cluster, is zero for p < Pc and strictly positive for p > Pc [1, 2]. Much 
is known rigorously and even more from computer simulations and renormalization group 
calculations, about the nature of the percolation transition in the independent case. In 
particular it is known rigorously that Pc is strictly greater than zero and less than one for 
d > 2 with Pc{d) a decreasing function of d, etc. We also know explicitly or have bounds 
for some of the various exponents associated with the divergence of different quantities, e.g. 
the mean finite cluster size, when p — > Pc- We even know exactly the scaling limit of the 
shape of the critical cluster on the triangular lattice [3] . 

It is generally believed that the critical properties, e.g. exponents, for independent per- 
colation, but not Pc, arc universal: they do not depend on the particular lattice but only on 
the dimensionality of the problem. The exponents arc also believed not to be changed when 
one considers systems for which the occupation probabilities for different sites are not inde- 
pendent, as long as the correlations between occupied sites decay rapidly, say exponentially 
[4]. 

Less is known about the percolation transition when there are long range correlations 
between occupied sites, e.g. when the correlations decay as a power law. Such power law 
decays occur in many physical systems and the nature of the percolation transition in such 
systems has come up recently in the study of two dimensional turbulence [5] , and of porous 
media, such as gels [6]. 

In a seminal work Weinrib and Halperin [4, 7] argued that the critical exponents of the 
percolation transition should depend only on the decay of the pair correlation G{r) in such 
systems. In particular for G{r) ~ r^" the transition should be in a universality class which 
depends only on a and d. Their analysis was based on considering the variance of the particle 
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density in a region of volume ^ , where ^ is the percolation correlation length which diverges 
as p y Pc- They found that if a < d these correlations are relevant if oi/ — 2 < 0. Here v is 
the critical exponent which describes the divergence of the percolation correlation length 
e.g. the average radius of gyration of the clusters in the independent percolation problem 
as p /" i-e. iip) ~ {pc — p)~^ ■ Weinrib and Halperin argued that systems that satisfy 
the above criteria belong to a new universality class for which the percolation correlation 
length exponent is viong — f [4, 7]. They also checked this using a renormalization group 
double expansion in e = 6 — d and in 5 = 4 — a. While the computations of WH were done 
only in the one loop approximation the exponent fiong was conjectured to be exact [4, 7] . 

As pointed out by WH their results are consistent with those based on renormalization 
group ideas, both in real and momentum space, on the percolation of like-pointing Ising 
type spins at the critical point, see [7] and references in there. There have also been some 
numerical tests of the WH predictions. For d — 2 Prakash et al. [8] have carried out Monte 
Carlo simulations for percolation on artificially generated power law correlated occupation 
probabilities on 1?. This study confirmed the predictions of Weinrib and Halperin. The 
only direct check of the WH prediction in d = 3 we are aware of is in [6] where the authors 
introduced a bond percolation model in 1?, called Pacman percolation. They argued that 
the pair correlation for their model decays as r~", with a close to 1, and obtained critical 
exponents which are consistent with WH, but since a was not known exactly the results are 
not fully conclusive. 

In this paper we study the percolation transition in three dimensions for two systems in 
which the long range correlations arise naturally from the microscopic dynamics: the mass- 
less harmonic crystal and the voter model on 1? . Both of these systems are known rigorously 
to have G{r) ~ r~^. They also have other similarities but are intrisicaly quite different. The 
existence and nature of the percolation transition in these systems is of interest in their own 
right. Using Monte Carlo simulations and finite size scaling we find the Pc for both models. 
We also find that both models have the same critical exponents as expected from the WH 
predictions of a long range percolation universality class. 

For the massless harmonic crystal in we define site x to be occupied if the scalar 
displacement field is greater than some prcassigncd value h and empty if < h. 
Percolation then corresponds to the existence of an infinite level set contour for (j){x) < h. 
The existence of percolation threshold, i.e. < Pc < 1; was proven by Bricmont, Lebowitz 
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and Maes [9] ior d — 3. There are however no previous calculations (known to us) concerning 
the actual value of Pc or of the critical exponents for this system. One expects intuitively 
that the Pc will be smaller than the Pc for independent percolation ,c.f. [8], but we know 
of no proof for this. Similarly a proof that Pc > for the harmonic crystal in (i > 3, or for 
the an-harmonic crystal in > 3 is still an open problem [10]. For d < 2, (j){x) is for any h, 
either plus or minus infinity, with probability 1, when the size of the system goes to infinity. 
Thus either all sites are occupied or all sites are empty. 

The voter model, often used for modeling various sociological and biological phenomena 
, is a lattice system in which a site x is occupied or empty according to whether the "voter" 
living there belongs to party A or B. Voters change their party affiliations according to a 
well defined stochastic dynamics [11]. The stationary state of this model is not known 
explicitly but many of its properties are known exactly. In particular it has many features 
in common with the harmonic crystal. Like the harmonic crystal, the stationary state of the 
voter model is trivial in d < 2; all sites occupied or all sites empty. On the other hand any 
p is possible on for d > 3, where the truncated pair correlation decays, as it does for the 
harmonic crystal, like r~^. No proof of the existence of a pc > is known for this system, 
i.e. the system could in principle percolate for arbitrary small p. For examples of systems 
where Pc < ^ for any e > see [12]. 

The outhne of the rest of the paper is as follows. In Section 2 we present the simulation 
methods and results the massless harmonic crystal. In particular we find Pc — 0.16 ± 0.01. 
In section 3 we study the voter model. We present a new efficient algorithm for simulat- 
ing this model and report the results from its implementation. We find in particular that 
Pc = 0.10 ± 0.01 compared with a — 0.16 obtained in [13] using a less reliable method. 
We conclude the paper with a brief discussion of some open problems. 

II. THE HARMONIC CRYSTAL 
A. Formulation 

Let X e designate the sites of a d-dimensional simple cubic lattice and (f){x) be the 
scalar displacement field at site x. The interaction potential in a box A with specified 
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boundary conditions (b.c), e.g. (f){x) — ior x on the boundary of A, has the form 

<x,y> x,y 

where J > and M > 0, < y > indicates nearest neighbor pairs, \x — y\ — 1, on TJ^. 
The sum is over all sites in A with the specified b.c. The Gibbs equihbrium distribution 
of the {0(x)} at a temperature = Z^p^ = e~^^ is then Gaussian with a 

covariance matrix f3C which is well defined for M > 0. 

The infinite volume limit Gibbs measure obtained when A is ,for M > 

,translation invariant, with < (f){x) >— and is independent of the boundary conditions 
[14]. When M ^ , //^ does not exist for d < 2 [14]. This is due to the fact that the 
fluctuations of the held, e.g < >, become unbounded for these dimensions. However, 

for d > 3 the Gibbs measure /j, obtained as the hmit of /^m when M ^ is well defined. (It 
is the same as the infinite volume limit of the measure in a box with M = and prescribed 
boundary values 0(a;) =0). In this limit the pair correlations between different sites have 
the long distance behavior r — \x — y\ for d > 2 [14]. 

Following [9] we define the occupation variable Ph{x) 

{1 a (f)(x) > h 
(2) 
ff (pix) < h 

and let p —< Ph{x) >^m, where the average is over the Gibbs measure /i^ . We can also 
define a new measure jj,^ on the occupation variables ph{x)— 0,1 by a projection of p^. All 
expectations involving a function of the occupation variables can be computed directly from 
/i^. The correlations between the occupation variables have the same asymptotic decay 
properties as those of the field variables 0, 

{Ph{x)ph{y))f,M - ~ I _ |rf_2 for d > 2 (3) 

\x y\ 

where ~ M"^ and the averages are with respect to p^{oT p^). In the limit M ^ 
the measure p, has a pair correlation that decays hke r^"*^ for d > 2. We note that p is not 
Gibbsian for any summable potential, c.f [15]. 
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B. Results 



Simulating the harmonic crystal on finite lattices is easy, the elements in a discrete Fourier 
transform of a harmonic crystal are independently distributed Gaussian random variables 
with easily computed variances [16]. We consider the system on a lattice with periodic 
boundary conditions and exclude the zero mode. This is essentially equivalent to fixing 
(0) = 0. 

There are many methods for obtaining the percolation threshold using data obtained 
from simulations on finite systems [1]. We used the method employed in [13, 17]. For a 
cube of linear size L let 




where rij is the number of clusters of j sites, defined by the occupation variables Phi^), and 
the average is taken over a large number of samples obtained from simulation of the model. 
We calculate for different sizes L and concentration of occupied sites p defined as in (2). 

One expects [1, 17, 18] that for large L,and (pc — p) <^ 1, Tl should have a finite size 
scaling form, 

L^'^Fl ~ L^F(L^(p — Pc)) + corrections to scaling, (5) 

where 7 is the critical exponent for the divergence as p y Pc of the second moment of the 
cluster size distribution, defined as the limit L ^ 00 of Corrections to scaling should 
go to zero for L — > 00. 

For p > Pc, for an infinite system the second moment of the cluster size distribution can 
be defined by excluding the infinite cluster. This diverges with a critical exponent 7' for 
p \ Pc- The finite system analog is ^ which is defined similarly to but not including 
the spanning cluster. scales as 

L~^r']^ ~ L^F\L~ [p — Pc)) + corrections to scaling. (6) 

It is believed that 7' = 7. 

According to finite size scaling theory the number of sites in the largest cluster in a finite 
system of linear size L, Pl{p), scales for |p — Pc| 1 as 

Pl{p) ~ L'^''^G{L^ {p — Pc)) + corrections to scaling (7) 
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[1, 18], where P is the critical exponent for the approach to zero of the fraction of sites 
belonging to the infinite cluster in an infinite system as p \ Pc- Using the hyper-scaling 
relation c? = 2^ + J we see that (5), (6) and (7) lead to the scahng form (5) being valid 
for all \p — Pc\ ^ 1 and large L. That is on a finite system we do not need to differentiate 
between p < pc or p > pc, we may include all the clusters when calculating Tl{p). 

Assuming (5) is valid for \p — pd <^ 1 the ratio Rl — ^ should become independent of 
L, for large L, at p — Pc- Plotting these ratios as a function of p for different sizes L and 
looking for the intersection of these different curves then yields Pc- The value of the ratios at 
the intersection point of the Rl curves should be equal to 2''+? giving us a way to measure 
^. Moreover ,we also have 




Thus the slopes of these curves should also give u. 

In Fig. 1 we present results of the simulation for the massless harmonic crystal on a 
cubic lattice with periodic boundary conditions. Each Tl was averaged over 48000 samples 
except for L — 160 where the average is over 2400 samples. To determine the error bars 
we have divided the output of the simulations into 10 parts and assuming that the averages 
are Gaussian distributed we evaluated the variance which we used as a measure of the 
uncertainty. From the intersection of the curves, after interpolation, we obtain Pc = 0.16 ± 
0.01. Comparing the slopes of the Rl curves for L = 80 and L = 40 we obtain u — 2.1± 0.5. 
From the value of Rl at the intersection point of the curves we obtain ^ = 1.8 ± 0.1. We 
actually computed for the sequences L = 10, 20, 40, 80, 160 and L = 15, 30, 60, 120. All 
the simulation results are consistent with what is plotted in Fig. 1 where we have used only 
part of these simulations since the plot is otherwise cluttered. These values clearly show 
that our system is in a different universality class from independent percolation since for the 
latter u = 0.876 ± 0.001 and ^ = 2.045 ± 0.001 [19]. 

The above method is good for finding the percolation threshold and the ratio of critical 
exponents ^ but clearly does not give good results for u. To obtain more precise result 
for the percolation correlation length exponent we evaluated the probability that there is a 
"wrapping cluster" , i.e. one that wraps around the torus, for different densities p of occupied 
sites and different linear sizes L. 

For fixed L we denote by p^-^^ the value of the density of occupied sites for which one half 
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of the realizations will have such a wrapping cluster. This should obey the following scaling 
relation ■p'^J^ — Pc ~ -^"^ [1]- For sizes between 30 and 100 we evaluated ■p^J^ from doing 
simulation in a range between p = 0.13 and p — 0.25 in steps of 0.005. For each such system 
24000 samples were generated. The slope of \o<g{jfJ^ — pc) versus log(L) should give us u. 
A plot of the results is presented in Fig. 2. The slope of the fitted straight line is 0.50 ± 0.01 
which gives u — 2.00 ±0.04. This is in good agreement with the theoretical prediction i/ — 2 
of Weinrib and Halperin [4] . 

We have used the obtained values of Pc, ^ and u to draw Fig. 3 where we see a good 
collapse of the data points to a smooth curve. 

We also calculated the ratio of the critical exponents ^. We did this by finding the 
fraction of sites that belong to the largest cluster in a system of linear size L, -^^|^^,when 
we simulate at the approximate critical density. From (7) we see that ^''^£2^^ ~ L~^. The 
result for systems of size from 40 to 170 averaged over 24000 samples is presented in (Fig. 4). 
From the slope of the fitted straight line we obtain ^ = 0.60 ± 0.01. Moreover, the fact that 
P{pc, L) follows well a power law behavior supports our contention that the true critical value 
is near pc = 0.16 ± 0.01. Observe also that 2^ + ^ = 3.0 ± 0.2 and thus the hyper-scaling 
relation is satisfied. 



III. THE VOTER MODEL 



Formulation 



Another system whose pair correlations decays like that of the massless harmonic crystal 
is the voter model in Z'' [11]. 

The voter model is defined through a stochastic time evolution. Each lattice site is 
occupied by a voter who can have two possible opinions, say yes or no. With rate the 
voter at site x adopts the opinion of one of his/her 2d neighbors chosen at random. More 
specifically letting p{x) — 0,1, x & Z'^, the time evolution of the voter model is specified by 
giving the rate C„(a:, p) for a change at site x when the configuration is given by p 



Cv{x,p) = - 

T 



1 - ^{2p{x) -1) J2 (My) - 1) 



\y-x\=l 
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where r sets the unit of time. 

It is clear that for the voter model on a finite set A C Z*^ with periodic or free boundary 
conditions (b.c), there will be only two possible stationary states: either p{x) = 1 or 
p{x) = for all X G A. The same is true for the voter model on an infinite lattice in one 
and two dimensions: the only stationary states are the consensus states. However for d > 3 
there are, as for the massless harmonic crystal, unique stationary states for every density p 
of positive spins, p — {p{x)). The correlations in this state decay as 

{p{x)p{y)) - = p(l - p)Gd{x - y) 

where Gd{x) is the probability for a random walker, starting ai x E Z'^ to hit the origin 
before escaping to infinity. It is well known that Gaix) ~ . i_2 for o? > 3, i.e the pair 

\x\ 

correlation for the voter model has the same long range behavior as the massless harmonic 
crystal. 

B. Simulation Method 

An efficient method to simulate the voter model is to consider a box Be of linear size C 
with stochastic boundary conditions, i.e. when a voter looks at the boundary he sees 1 with 
probability p and with probability 1 — p. It is then possible to show that the distribution 
of the configuration of voters in a box of size L < C centered inside Be and far away 
from the boundary will approach the steady state measure (restricted to Bl) with density p 
for the voter model when C ^ oo. In order to sample from the measure for the voter model 
inside Be with such stochastic boundary conditions we use the following algorithm: Start a 
random walk from each site oi B^ and let these random walks move independently until two 
of them meet in which case they coalesce. When a random walk hits the boundary of Be it 
is frozen. We continue this until all the random walkers either coalesce or are frozen. After 
this is done we independently for each frozen walker, assign the value 1 with probability p 
and the value with probability 1 — p, then assign that same value to its ancestors, that is 
all the random walkers that have coalesced with it. In this way we assign values 1 or zero to 
all the sites in i^^. One can prove that in this way we sample configurations inside B^ with 
the distribution coming from the voter model in Be with the stochastic boundary conditions 
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discussed above. The advantage of this way of simulating is that one is guaranteed that the 
samphng is from the steady state measure with these boundary conditions. 

C. Results 

Using this method of generating configurations inside Bl for different p we looked for 
a spanning cluster inside Bl- We did simulations for sizes L — 10, 15, 20, 25 and 30 with 
jC — 160. The results which are the same for all C in the range (120, 160) are presented in 
Fig. 5. If we assume the scaling form for the spanning probability [1] 

Ili^{p)^F{{p-p,)Lyn (9) 

then by collapsing the data Fig. 6 we obtain = 0.10 ± 0.01 and z/ = 2 ± 0.2. 

To find ^ we measured ^ and we assume the scaling form (5). Note that in this case 
we do not have periodic boundary conditions. Results from the simulation are presented in 
Fig. 7. Collapsing the data Fig. 8 we obtain p^ = 0.10 ± 0.01, ^ = 1.9 ± 0.2 and i/ = 2 ± 0.2. 

Analogous simulation measurements for P{p,L) gave ^ = 0.6 ± 0.1. As in the case of 
the massless harmonic crystal the exponents we found satisfy the hypcr-scaling relation 
2^ + ^ = d. The exponents for both the massless harmonic crystal and voter model seem 
to agree within the error bars. 

D. Comparison of Pc with previous simulations 

The percolation transition in the d — 3 voter model was first investigated in [13]. This 
was done by considering voters who occasionally change their opinions spontaneously, i.e. 
independently of what their neighbors are doing. They do this with probability A. In terms 
of flip rates one has 

C{x,p) = (1 - X)C,{x,p) + ^ [1 + (1 - 2p){2p{x) - 1)] , 

T 

where < p < 1 and < A < 1 and C„ is the voter model flip rates. This leads to a sta- 
tionary state in any periodic box of size L'^ with density of pluses equal to p. As A increases 
from to 1 we go from the voter model to an independent flip model. The stationary state 
of the latter is a product measure with density p. This model was studied rigorously in [20] 
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where it was named the noisy voter model. 

In [13] the authors used (5), on simulation results of the noisy voter model on lattices 
with periodic boundary conditions, to obtain Pc{X) for A > 0.1. For d — 3 they found by 
extrapolation Pc(^) ~ 0.16 as A ^ 0. 

We have repeated the simulations in [13] for larger lattice sizes and smaller values of A. 
We simulated systems with A as small as 0.01 each with 24000 "effectively uncorrelated" 
samples and sizes up to 80. Prom our results we can extrapolate Pc(A) 0.15 as A — > 0, a 
value slightly lower than the result in [13]. We also observed that, as expected, the critical 
exponents for the noisy voter model agree, for the given range of A, with the critical expo- 
nents of independent percolation. 

This leaves a significant difference with the result for Pc obtained in the previous section. 
We believe that the answer lies in the necessary extrapolation to A = 0. Since the autocorre- 
lation time grows exponentially with lambda, this means we have to wait for more and more 
Monte Carlo steps to get independent samples. To check this explanation we investigated 
the percolation transition in the harmonic crystal with a mass M. This mass acts much like 
the random flips in the voter model. For both models the pair correlation decays exponen- 
tially. In the harmonic crystal the characteristic length scale is C,m = ji- An easy calculation 



shows that the characteristic length scale for the noisy voter model is = y The noisy 
voter model with the smallest lambda that we simulated , A = 0.01, thus corresponds to 
roughly equal to 4 (unit distance is the lattice spacing) . In the language of the massive 
harmonic crystal this corresponds to M ~ 0.25. Estimating the percolation threshold of the 
massless harmonic crystal by the extrapolation method we used for the voter model using 
M > 0.25 yields a Pc{M) ~ 0.21 when M — 0. This is obviously a large overestimate of 
Pc = 0.16 which was obtained by directly simulating the massless harmonic crystal. This 
shows that the extrapolation method greatly overestimates the true Pc- 

IV. CONCLUDING REMARKS 

We have performed Monte Carlo simulations to obtain the critical percolation density and 

some critical exponents for the massless harmonic crystal and the voter model in Z^. We 
found, for the first time a value of Pc for the former and using a novel method of simulation 
for the voter model found a new more reliable value of Pc for this system. The critical 
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exponents for both models agree within the error bars. This suggests that both percolation 
models are in the same universality class and confirms the theoretical predictions made in 
[4]. The result for the correlation length critical exponent v — 2 supports the conjecture by 

WH that the relation v = - is exact. 

a 

It is believed that not only the critical exponents but also the finite size scaling functions 
are universal. While this is certainly consistent with our simulations we have not checked 
this carefully. Such a check would require measuring quantities for the two systems in the 
same way. This is not what we have done here as we wanted to use the most efficient method 
for each system. 

We mention here that there has been much activity in generalizing the voter model in 
various ways [21]. Based on our present work we expect that the nature of the percolation 
transition in these models will depend only on the asymptotic behavior of G{r). We have 
however not investigated this. Our simulation method may be extendable to some of these 
systems. 

The reported simulations were done on a Sun Microsystems HPC- 10000 system. 
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FIG. 3: Plot of riL~"~^ versus {p — Pc)L'i^ for the d = 3 massless harmonic crystal for = 0.16, 
1/ = 2 and ^ = 1.8. We have plotted data points for L = 30, 60, 120 for p = 0.13 to 0.16 in steps of 
0.005 and for L = 40, 80, 160 for p = 0.13 to 0.18 in steps of 0.005. 
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FIG. 4: Plot of log(P(pc, L)) versus log(L). The slope of the straight line gives ^ = 0.60 ± 0.01 
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FIG. 5: Plot of Hi versus p for the d = 3 voter model. We have plotted data points for L = 
10, 15, 20, 25 and 30 from p = 0.06 to p = 0.205 in steps of 0.005. Each point is an average over 
10^ samples. The error bars are not shown since on this scale they are too small. 
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FIG. 6: Plot of n^, versus {p — pc)Lv for the d = 3 voter model for pc = 0.105 and u = 2. We have 
used the same data that was used to create Fig. 5 
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FIG. 7: Plot of TlL~'^ versus p for the d = 3 voter model. We have plotted data points for 

L = 10, 15, 20, 25 and 30 from p = 0.06 to p = 0.205 in steps of 0.005. Each point is an average 
over 10^ samples. The error bars are not shown since on this scale they are too small. 
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FIG. 8: Plot of FlL 1 versus {p - Pc)lI for the d = 3 voter model for pc = 0.105, v = 2 and 
^ = 1.9. We have used the same data that was used to create Fig. 7 
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